data {
  
  // number of total observations -- i.e. all conjoint profiles that were seen
  // if there are n_resp respondents who each view
  // n_plans plans, then n = n_resp * n_plans
  int<lower=1> n;
  
  // binary vector of length n indicating whether a plan was preferred
  real<lower=0,upper=1> y[n]; 

  // table linking y (outcomes) to respondents, for random effects
  // n_resp = number of respondents in the survey
  // ii[j] gives the respondent index for observation j. ii has the same
  // length as the response vector y.
  int<lower=1> n_resp;
  int<lower=1,upper=n_resp> ii[n];

  // levels of conjoint
  // k_factor = number of levels for a given factor, excluding baseline level
  // factor_mat = binary matrices indicating which levels of a factor were
  // "turned on" for a given observation. 
  int<lower=1> k_educ;
  matrix[n,k_educ] educ_mat;

  int<lower=1> k_hieduc;
  matrix[n,k_hieduc] hieduc_mat;

  int<lower=1> k_invest;
  matrix[n,k_invest] invest_mat;

  int<lower=1> k_gov;
  matrix[n,k_gov] gov_mat;

  int<lower=1> k_workers;
  matrix[n,k_workers] workers_mat;

  int<lower=1> k_local;
  matrix[n,k_local] local_mat;
  
  // individual-level covariates
  // k_indvl = number of individual-level covariates
  // ind_covariates = matrix of individual-level covs
  int<lower=1> k_indvl; 
  matrix[n_resp,k_indvl] ind_covariates;
}

parameters {
  
  // beta parameter vectors - one for each respondent. these will be
  // given N(0, 1) priors, so they need to be transformed to have 
  // the correct SD's before entering the likelihood. See the 
  // transformed parameters block below. 
  matrix[n_resp, k_educ] beta_educ_tilde;
  matrix[n_resp, k_hieduc] beta_hieduc_tilde;
  matrix[n_resp, k_invest] beta_invest_tilde;
  matrix[n_resp, k_gov] beta_gov_tilde;
  matrix[n_resp, k_workers] beta_workers_tilde;
  matrix[n_resp, k_local] beta_local_tilde;

  
  // gamma params for hierarchical mean of beta params
  matrix[k_educ, k_indvl] gamma_educ;
  matrix[k_hieduc, k_indvl] gamma_hieduc;
  matrix[k_invest, k_indvl] gamma_invest;
  matrix[k_gov, k_indvl] gamma_gov;
  matrix[k_workers, k_indvl] gamma_workers;
  matrix[k_local, k_indvl] gamma_local;

  // intercept term (alpha_i in paper's notation)
  vector[n_resp] intercept_tilde; 

  // second-level parameters for hierarchical mean of intercept
  row_vector[k_indvl] gamma_intercept;
  

  // Standard deviation of random effects
  real<lower=0> intercept_sd;
  real<lower=0> educ_sd[k_educ];
  real<lower=0> hieduc_sd[k_hieduc];
  real<lower=0> invest_sd[k_invest];
  real<lower=0> gov_sd[k_gov];
  real<lower=0> workers_sd[k_workers];
  real<lower=0> local_sd[k_local];


  // Residual standard deviation
  real<lower=0> residual_sd; 

}


transformed parameters{
  matrix[n_resp, k_educ] beta_educ;
  matrix[n_resp, k_hieduc] beta_hieduc;
  matrix[n_resp, k_invest] beta_invest;
  matrix[n_resp, k_gov] beta_gov;
  matrix[n_resp, k_workers] beta_workers;
  matrix[n_resp, k_local] beta_local;


  vector[n_resp] intercept; 

  // define random effects that enter data likelihood
  for (i in 1:n_resp){
    for (j in 1:k_educ){
      beta_educ[i, j] =  dot_product(ind_covariates[i, ], gamma_educ[j, ]) + educ_sd[j] * beta_educ_tilde[i,j]; 
    }
  }

  for (i in 1:n_resp){
    for (j in 1:k_hieduc){
      beta_hieduc[i, j] =  dot_product(ind_covariates[i, ], gamma_hieduc[j, ]) + hieduc_sd[j] * beta_hieduc_tilde[i,j]; 
    }
  }  
  
  for (i in 1:n_resp){
    for (j in 1:k_invest){
      beta_invest[i, j] =  dot_product(ind_covariates[i, ], gamma_invest[j, ]) + invest_sd[j] * beta_invest_tilde[i,j]; 
    }
  }
  for (i in 1:n_resp){
    for (j in 1:k_gov){
      beta_gov[i, j] =  dot_product(ind_covariates[i, ], gamma_gov[j, ]) + gov_sd[j] * beta_gov_tilde[i,j]; 
    }
  }
  for (i in 1:n_resp){
    for (j in 1:k_workers){
      beta_workers[i, j] =  dot_product(ind_covariates[i, ], gamma_workers[j, ]) + workers_sd[j] * beta_workers_tilde[i,j]; 
    }
  }
  for (i in 1:n_resp){
    for (j in 1:k_local){
      beta_local[i, j] =  dot_product(ind_covariates[i, ], gamma_local[j, ]) + local_sd[j] * beta_local_tilde[i,j]; 
    }
  }
  for (i in 1:n_resp){
    intercept[i] = dot_product(gamma_intercept, ind_covariates[i,]) + intercept_sd * intercept_tilde[i];
  }
}


// priors and likelihood 
model {

  vector[n] lin_pred;

  

  // prior on (transformed) random effects 
  for (i in 1:n_resp){
    for (j in 1:k_educ){
        target += normal_lpdf(beta_educ_tilde[i, j] | 0, 1);    
     }
    for (j in 1:k_hieduc){
        target += normal_lpdf(beta_hieduc_tilde[i, j] | 0, 1);    
     }
    for (j in 1:k_invest){
        target += normal_lpdf(beta_invest_tilde[i, j] | 0, 1);    
     }
    for (j in 1:k_gov){
        target += normal_lpdf(beta_gov_tilde[i, j] | 0, 1);    
     }
    for (j in 1:k_workers){
        target += normal_lpdf(beta_workers_tilde[i, j] | 0, 1);    
     }
    for (j in 1:k_local){
        target += normal_lpdf(beta_local_tilde[i, j] | 0, 1);    
     }
     target += normal_lpdf(intercept_tilde[i] | 0, 1);
  }

  

  // priors on hierarchical coefficients and standard deviations
  for (i in 1:k_educ){
    gamma_educ[i] ~ normal(0, 10);
    educ_sd[i] ~ cauchy(0, 2);
  }
  for (i in 1:k_hieduc){
    gamma_hieduc[i] ~ normal(0, 10);
    hieduc_sd[i] ~ cauchy(0, 2);
  } 
  for (i in 1:k_invest){
    gamma_invest[i] ~ normal(0, 10);
    invest_sd[i] ~ cauchy(0, 2);
  }
  for (i in 1:k_gov){
    gamma_gov[i] ~ normal(0, 10);
    gov_sd[i] ~ cauchy(0, 2);
  }
  for (i in 1:k_workers){
    gamma_workers[i] ~ normal(0, 10);
    workers_sd[i] ~ cauchy(0, 2);
  }
  for (i in 1:k_local){
    gamma_local[i] ~ normal(0, 10);
    local_sd[i] ~ cauchy(0, 2);
  }

  // intercept standard deviation prior 
  intercept_sd ~ cauchy(0, 2);

  // residual sd prior
  residual_sd ~ cauchy(0, 2);


 

  // DATA MODEL // 

  // get linear predictor 
  for (i in 1:n) {
    lin_pred[i] = dot_product(beta_educ[ii[i], ],  educ_mat[i, ]) + 
            dot_product(beta_hieduc[ii[i], ],  hieduc_mat[i, ]) + 
            dot_product(beta_invest[ii[i], ],  invest_mat[i, ]) + 
            dot_product(beta_gov[ii[i], ],  gov_mat[i, ]) + 
            dot_product(beta_workers[ii[i], ],  workers_mat[i, ]) + 
            dot_product(beta_local[ii[i], ],  local_mat[i, ]) + 
            intercept[ii[i]];
  }

  // sampling statement for the outcome variables y ~ N(\beta_i * X_i, residual_sd)
  target += normal_lpdf(y | lin_pred, residual_sd);
}


